home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 2.iso / STUTTGART / LANG / HUGS1 / hs / matrix < prev    next >
Text File  |  1995-02-14  |  3KB  |  94 lines

  1. -- Some simple Hugs programs for manipulating matrices.
  2. --
  3.  
  4. type Matrix k = [Row k]          -- matrix represented by a list of its rows
  5. type Row k    = [k]              -- a row represented by a list of literals
  6.  
  7. -- General utility functions:
  8.  
  9. shapeMat    :: Matrix k -> (Int, Int)
  10. shapeMat mat = (rows mat, cols mat)
  11.  
  12. rows        :: Matrix k -> Int
  13. rows mat     = length mat
  14.  
  15. cols        :: Matrix k -> Int
  16. cols mat     = length (head mat)
  17.  
  18. idMat       :: Int -> Matrix Int
  19. idMat 0      = []
  20. idMat (n+1)  = [1:copy n 0] ++ map (0:) (idMat n)
  21.  
  22. -- Matrix multiplication:
  23.  
  24. multiplyMat                     :: Matrix Int -> Matrix Int -> Matrix Int
  25. multiplyMat a b | cols a==rows b = [[row `dot` col | col<-b'] | row<-a]
  26.                 | otherwise      = error "incompatible matrices"
  27.                  where v `dot` w = sum (zipWith (*) v w)
  28.                        b'        = transpose b
  29.  
  30. -- An attempt to implement the standard algorithm for converting a matrix
  31. -- to echelon form...
  32.  
  33. echelon   :: Matrix Int -> Matrix Int
  34. echelon rs
  35.     | null rs || null (head rs) = rs
  36.     | null rs2                  = map (0:) (echelon (map tail rs))
  37.     | otherwise                 = piv : map (0:) (echelon rs')
  38.       where rs'            = map (adjust piv) (rs1++rs3)
  39.             (rs1,rs2)      = span leadZero rs
  40.             leadZero (n:_) = n==0
  41.             (piv:rs3)      = rs2
  42.  
  43. -- To find the echelon form of a matrix represented by a list of rows rs:
  44. -- 
  45. -- {first line in definition of echelon}:
  46. --  If either the number of rows or the number of columns in the matrix
  47. --  is zero (i.e. if null rs || null (head rs)), then the matrix is
  48. --  already in echelon form.
  49. -- 
  50. -- {definition of rs1, rs2, leadZero in where clause}:
  51. --  Otherwise, split the matrix into two submatrices rs1 and rs2 such that
  52. --  rs1 ++ rs2 == rs  and all of the rows in rs1 begin with a zero.
  53. --
  54. -- {second line in definition of echelon}:
  55. --  If rs2 is empty (i.e. if null rs2) then every row begins with a zero
  56. --  and the echelon form of rs can be found by adding a zero on to the
  57. --  front of each row in the echelon form of (map tail rs).
  58. --
  59. -- {Third line in definition of echelon, and definition of piv, rs3}:
  60. --  Otherwise, the first row of rs2 (denoted piv) contains a non-zero
  61. --  leading coefficient.  After moving this row to the top of the matrix
  62. --  the original matrix becomes  piv:(rs1++rs3).
  63. --  By subtracting suitable multiples of piv from (suitable multiples of)
  64. --  each row in (rs1++rs3) {see definition of adjust below}, we obtain a
  65. --  matrix of the form:
  66. --
  67. --          <----- piv ------>
  68. --          __________________
  69. --          0  |
  70. --          .  |
  71. --          .  |      rs'        where rs' = map (adjust piv) (rs1++rs3)
  72. --          .  |
  73. --          0  |
  74. --
  75. --  whose echelon form is  piv : map (0:) (echelon rs').
  76. --
  77.  
  78. adjust              :: Num a => Row a -> Row a -> Row a
  79. adjust (m:ms) (n:ns) = zipWith (-) (map (n*) ms) (map (m*) ns)
  80.  
  81. -- A more specialised version of this, for matrices of integers, uses the
  82. -- greatest common divisor function gcd in an attempt to try and avoid
  83. -- result matrices with very large coefficients:
  84. --
  85. -- (I'm not sure this is really worth the trouble!)
  86.  
  87. adjust'              :: Row Int -> Row Int -> Row Int
  88. adjust' (m:ms) (n:ns) = if g==0 then ns
  89.                                 else zipWith (\x y -> b*y - a*x) ms ns
  90.                         where g = gcd m n
  91.                               a = n `div` g
  92.                               b = m `div` g
  93. -- end!!
  94.